## This program generates the top panels of Figure 5 from 
##  "Financial Constraints, Sectoral Heterogeneity, and the Cyclicality of Investment" by Cooper Howes


##### Step 0: Preliminaries #####


## Clear variables and data
rm(list=ls())

## Load packages
library(sandwich)
library(tidyverse)
library(dynlm)
library(glue)


##### Step 1: Load the data #####


## BEA fixed asset data by industry are downloaded here: https://www.bea.gov/itable/fixed-assets
##  The .csv files used in the published version are in the "BEA fixed asset files" folder

#--- Real investment ---#

## Table 3.8ESI
real_fa <- t(read.csv("BEA fixed asset files/BEA_FA_invest.csv",skip=4)[1:78,-(1:2)])

## Table 3.8E
real_equip <- t(read.csv("BEA fixed asset files/BEA_equip_invest.csv",skip=4)[1:78,-(1:2)])

## Table 3.8S
real_struct <- t(read.csv("BEA fixed asset files/BEA_structures_invest.csv",skip=4)[1:78,-(1:2)])

## Table 3.8I
real_ip <- t(read.csv("BEA fixed asset files/BEA_IP_invest.csv",skip=4)[1:78,-(1:2)])


#--- Nominal investment ---#

## Table 3.7ESI
nominal_fa <- t(read.csv("BEA fixed asset files/BEA_FA_invest_nominal.csv",skip=4)[1:78,-(1:2)])

## Table 3.7E
nominal_equip <- t(read.csv("BEA fixed asset files/BEA_equip_invest_nominal.csv",skip=4)[1:78,-(1:2)])

## Table 3.7S
nominal_struct <- t(read.csv("BEA fixed asset files/BEA_structures_invest_nominal.csv",skip=4)[1:78,-(1:2)])

## Table 3.7I
nominal_ip <- t(read.csv("BEA fixed asset files/BEA_ip_invest_nominal.csv",skip=4)[1:78,-(1:2)])


## Create price indices by dividing nominals by reals
prices_fa <- nominal_fa/real_fa*100
prices_equip <- nominal_equip/real_equip*100
prices_struct <- nominal_struct/real_struct*100
prices_ip <- nominal_ip/real_ip*100

## Extract industry names
industries <- as.character(read.csv("BEA fixed asset files/BEA_FA_invest.csv",skip=4)[1:78,2])

## Import monetary policy shocks
mp_shocks <- ts(read.csv("Aggregate_data.csv"),
                start=c(1966,1),end=c(2018,1),frequency=4)[,"rr_shock"]

macronames <- colnames(mp_shocks)

## Pick subset of industries for which to run regressions 
## Note: this defaults to including several industries not shown in the paper
istar <- c(1,5,9,10,11,12,24,33,34,35,44,49,55,66,67,72,75,58)

## Extract names to be plotted
outcomenames <- c("rr_shock",industries)
plotters <- industries[istar]


##### Step 2: Convert data from annual to quarterly #####


## The BEA fixed asset series will be merged with the MP shock series, which are
##  quarterly. To combine them, I first turn the annual series into quarterly by
##  filling in all quarters with the annual observation. Then, before running the
##  regressions, I aggregate everything back to annual frequency by taking averages.

## Create empty matrix
zeromatrix <- matrix(0,nrow=dim(real_fa)[1]*4,ncol=dim(real_fa)[2])

fa_ts <- zeromatrix
equip_ts <- zeromatrix
struct_ts <- zeromatrix
ip_ts <- zeromatrix

fa_price_ts <- zeromatrix
equip_price_ts <- zeromatrix
struct_price_ts <- zeromatrix
ip_price_ts <- zeromatrix

fa_nom_ts <- zeromatrix
equip_nom_ts <- zeromatrix
struct_nom_ts <- zeromatrix
ip_nom_ts <- zeromatrix


## Fill in annual data for each quarter
for(i in 1:(dim(real_fa)[1])){
  fa_ts[(4*(i-1)+1):(4*i),] <-matrix(rep(real_fa[i,],4),nrow=4,byrow=T)
  equip_ts[(4*(i-1)+1):(4*i),] <-matrix(rep(real_equip[i,],4),nrow=4,byrow=T)
  struct_ts[(4*(i-1)+1):(4*i),] <-matrix(rep(real_struct[i,],4),nrow=4,byrow=T)
  ip_ts[(4*(i-1)+1):(4*i),] <-matrix(rep(real_ip[i,],4),nrow=4,byrow=T)
  
  fa_price_ts[(4*(i-1)+1):(4*i),] <-matrix(rep(prices_fa[i,],4),nrow=4,byrow=T)
  equip_price_ts[(4*(i-1)+1):(4*i),] <-matrix(rep(prices_equip[i,],4),nrow=4,byrow=T)
  struct_price_ts[(4*(i-1)+1):(4*i),] <-matrix(rep(prices_struct[i,],4),nrow=4,byrow=T)
  ip_price_ts[(4*(i-1)+1):(4*i),] <-matrix(rep(prices_ip[i,],4),nrow=4,byrow=T)
  
  fa_nom_ts[(4*(i-1)+1):(4*i),] <-matrix(rep(nominal_fa[i,],4),nrow=4,byrow=T)
  equip_nom_ts[(4*(i-1)+1):(4*i),] <-matrix(rep(nominal_equip[i,],4),nrow=4,byrow=T)
  struct_nom_ts[(4*(i-1)+1):(4*i),] <-matrix(rep(nominal_struct[i,],4),nrow=4,byrow=T)
  ip_nom_ts[(4*(i-1)+1):(4*i),] <-matrix(rep(nominal_ip[i,],4),nrow=4,byrow=T)
}


## Convert to time series
fa_ts <- ts(fa_ts,start=c(1947,1),freq=4)
struct_ts <- ts(struct_ts,start=c(1947,1),freq=4)
equip_ts <- ts(equip_ts,start=c(1947,1),freq=4)
ip_ts <- ts(ip_ts,start=c(1947,1),freq=4)

fa_price_ts <- ts(fa_price_ts,start=c(1947,1),freq=4)
struct_price_ts <- ts(struct_price_ts,start=c(1947,1),freq=4)
equip_price_ts <- ts(equip_price_ts,start=c(1947,1),freq=4)
ip_price_ts <- ts(equip_price_ts,start=c(1947,1),freq=4)

fa_nom_ts <- ts(fa_nom_ts,start=c(1947,1),freq=4)
struct_nom_ts <- ts(struct_nom_ts,start=c(1947,1),freq=4)
equip_nom_ts <- ts(equip_nom_ts,start=c(1947,1),freq=4)
ip_nom_ts <- ts(equip_nom_ts,start=c(1947,1),freq=4)


## Combine MP shocks with fixed asset series
fa_data <- ts.intersect(mp_shocks,fa_ts)
equip_data <- ts.intersect(mp_shocks,equip_ts)
struct_data <- ts.intersect(mp_shocks,struct_ts)
ip_data <- ts.intersect(mp_shocks,ip_ts)

fa_price_data <- ts.intersect(mp_shocks,fa_price_ts)
equip_price_data <- ts.intersect(mp_shocks,equip_price_ts)
struct_price_data <- ts.intersect(mp_shocks,struct_price_ts)
ip_price_data <- ts.intersect(mp_shocks,ip_price_ts)

fa_nom_data <- ts.intersect(mp_shocks,fa_nom_ts)
equip_nom_data <- ts.intersect(mp_shocks,equip_nom_ts)
struct_nom_data <- ts.intersect(mp_shocks,struct_nom_ts)
ip_nom_data <- ts.intersect(mp_shocks,ip_nom_ts)




##### Step 3: Run regressions and store coefficients #####



## Set length of response horizon (in years)
numres <- 4

## Set number of standard deviations to show in errors (default: 1.65)
se_scale <- 1.65

## Set number of autoregressive lags
autolags <- 4

## Set beginning and ending dates for regressions
begin <- c(1970)
fin   <- c(2008)

## Specify the type of MP shock to be used in regressions (default: "rr")
shocktype <- "rr"

## Specify regression formula
f <- formula( glue("L(y,-i) ~ trend(y) + {shocktype}_shock +L(y,1:autolags) ")  ) 

## Create empty arrays to hold coefficients and standard errors
gamma    <- matrix(0,nrow=(numres),ncol=length(outcomenames))
gamma_se <- matrix(0,nrow=(numres),ncol=length(outcomenames))

fa_gamma    <- matrix(0,nrow=(numres),ncol=length(outcomenames))
fa_gamma_se <- matrix(0,nrow=(numres),ncol=length(outcomenames))

equip_gamma    <- matrix(0,nrow=(numres),ncol=length(outcomenames))
equip_gamma_se <- matrix(0,nrow=(numres),ncol=length(outcomenames))

struct_gamma    <- matrix(0,nrow=(numres),ncol=length(outcomenames))
struct_gamma_se <- matrix(0,nrow=(numres),ncol=length(outcomenames))

ip_gamma    <- matrix(0,nrow=(numres),ncol=length(outcomenames))
ip_gamma_se <- matrix(0,nrow=(numres),ncol=length(outcomenames))


## Outer loop over each type of fixed assets (total, equipment, structures, or intellectual property)
for(asset_type in c("fa","equip","struct","ip")){
  
  eval(parse(text=paste0("data <- ",asset_type,"_data")))

  colnames(data) <- outcomenames
  plotvars <- match(plotters,colnames(data))
  
  ## Inner loop over industries
  for(j in plotvars){
    
    n <- outcomenames[j]
    
    y <- log(data[,n])
    
    
    ## Innermost loop over response horizons
    for(i in c(1:numres)){
      
      ## Collapse data to annual frequency before running regressions
      data %>% aggregate.ts(.,frequency=4,FUN=sum)
      
      ## Set last year so that we don't use any data post-2008
      endyr <- 2008-i
      
      ## Estimate regression and calculate Newey-West standard errors
      eq <- dynlm(formula=f,start=begin,end=endyr,data=data)
      newey <- NeweyWest(eq,prewhite=FALSE)
      
      ## Save coefficient estimates and standard errors
      gamma[i,j] <- coef(summary(eq))[paste0(shocktype,"_shock"),"Estimate"]
      gamma_se[i,j] <- sqrt(newey[paste0(shocktype,"_shock"),paste0(shocktype,"_shock")])
      
      
    }
    
  }
  
  ## Save coefficient estates
  eval(parse(text=paste0(asset_type,"_gamma <- gamma")))
  eval(parse(text=paste0(asset_type,"_gamma_se <- gamma_se")))
  
}




##### Step 4: Plot impulse responses shown in top panel of Figure 5 #####



## Create response horizon over which coefficients will be plotted
years <- 1:4

## Full names for types of fixed assets to be used in plots
longtypenames <- c("Total fixed asset investment","Equipment investment","Structures investment")

## Create new window
windows(width=12,height=4)
par(oma=c(0,0,0,0))
m <- matrix(1:3,nrow=1,byrow = TRUE)
layout(mat = m)

## Create plot index to loop over
plotindex <- 1

## Plot figures for reach type of fixed asset
for(types in c("fa","equip","struct")){
  
  ## Replace "gamma" with the relevant set of coefficients to be plotted
  eval(parse(text=paste0("gamma <- ",types,"_gamma")))
  eval(parse(text=paste0("gamma_se <- ",types,"_gamma_se")))
  
  ## Specify which industries (indexed by their position in istar) to show
  inumbers <- c(5,4,9,10,11,13,14,18,1)
  
  ## Define which variables to plot (excluding the total, which is plotted separately)
  multiplots_all <- match(plotters[inumbers],colnames(data))
  multiplots <- multiplots_all[-c(1,length(inumbers))]
  
  ## Choose colors to use for plots
  plotcolors <- c("red","forestgreen","orange","purple","salmon",
                  "green2","deeppink2")
  
  ## Specify y-axis range
  yrange <- c(-0.05,0.05)
  
  ## Placeholder to loop through list of colors
  tempcolor <- 1
  
  ## Create plot with first industry
  par(mar=c(4,5,2,1))
  plot(x=years,y=gamma[,multiplots_all[1]],type="b",col="blue",lwd=4,pch=1,ylim=yrange,
       axes=F,frame=T,ann=FALSE,cex.main=1.25,cex=2)
  axis(1,cex.axis=2,at=years,lab=years)
  axis(2,cex.axis=2,at=pretty(yrange),lab=paste0(pretty(yrange)*100, "%"),las=TRUE)
  if(types=="equip"){mtext("Years after shock",side=1,line=2.5,cex=1.25)}
  abline(a=0,b=0,col="black",lwd=1)
  
  ## Loop to add lines for other industries
  for(i in multiplots){
    
    lines(x=years,y=gamma[,i],type="b",col=plotcolors[tempcolor],lwd=4,pch=(tempcolor+1),cex=2)
    
    ## Move to next color
    tempcolor <- tempcolor+1
  }
  
  ## Add line for total economy and title
  lines(x=years,y=gamma[,plotvars[1]],type="l",col="black",lwd=6)
  title(main=longtypenames[plotindex], col.main="blue", font.main=2,cex.main=1.75)
  
  ## Move to next index for the list of industries to plot
  plotindex <- plotindex+1
}


